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abstract 


A  computational  method  of  free-surface  wave  height  prediction 
for  a  ship  in  steady  motion  has  been  developed  through  application 
of  Michell's  thin  ship  theory.  Numerical  computation  for  two 
different  hulls  has  been  presented.  The  results  for  the  Wigley 
hull  agree  well  with  the  experiment.  However,  the  agreement  between 
the  numerical  results  and  the  experiment  for  the  SL7  container  ship 
is  not  satisfactory.  Further  study  in  analytical  method  is  necessary 
to  improve  the  numerical  results  for  a  ship  with  a  bulbous  bow. 


ADMINISTRATIVE  INFORMATION 

This  work  was  performed  under  the  General  Hydromechanics  Research  Program  and 
was  authorized  by  the  Naval  Sea  Systems  Command  (NAVSEA) ,  Hull  Research  and  Tech¬ 
nology  Office.  Funding  was  provided  under  Program  Element  61153N,  Task  Area 
SR  0230101,  and  Work  Unit  1542-700. 


INTRODUCTION 

The  analytical  method  for  computing  the  free-surface  wave  height  of  a  ship  is 
one  of  the  most  important  tasks  in  naval  hydrodynamics.  The  steady  ship  wave  pro¬ 
file  is  usually  linearly  superimposed  on  other  free-surface  wave  profiles.  In  the 
computation  of  the  relative  bow  motion  of  a  naval  ship  in  waves,  it  is  necessary  to 
know  the  wave  height  due  to  the  steady  forward  motion  before  computing  the  transfer 
functions.  The  analytical  study  of  deck  wetness  also  requires  the  knowledge  of  the 
wave  profiles  in  steady  motion. 

The  present  computational  method  has  been  developed  through  an  application  of 
Michell's  thin  ship  theory.  This  theory  is  first-order  and  linearized.  Even  though 
there  is  some  limitation  in  the  application  of  thin  ship  theory  to  ships  with  large 
block  coefficients  and  with  blunt  bows,  the  results  of  this  theory  are  reliable  for 
various  hull  forms  as  concluded  in  Referetrtee  1.  .  -  *  ' 

Two  hull  forms  have  been  chosen  for  the  numerical  computations  in  this  presen¬ 
tation:  a  Wigley  hull  and  an  SL7  container  ship.  The  numerical  results  for  the 

Wigley  hull,  which  is  a  mathematical  hull  form,  show  satisfactory  agreement  with  the 
experiment.  However,  the  numerical  results  for  the  SL7  are  not  satisfactory.  There 
are  some  discrepancies  near  the  bow  region  for  even  small  Froude  numbers. 


*A  complete  listing  of  references  is  given  on  page  17. 
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EQUATIONS  OF  FREE-SURFACE  WAVE  HEIGHT 

The  coordinate  system  o,x,y,z  is  fixed  in  the  ship  with  o,x,y  containing  the 

longitudinal  centerplane,  o,y,z  the  midship  section,  and  o,x,z  the  undisturbed 

water  surface  if  the  ship  is  at  rest.  The  vector  o,x  points  toward  the  stern,  and 

the  ship  is  moving  with  speed  U  toward  the  negative  x-axis.  With  assumptions  that 

the  fluid  is  inviscid  and  the  flow  is  irrotational,  the  governing  equation  and 

2 

boundary  conditions  are: 

(i)  in  the  fluid  domain 


<b  +  <b  +  d>  =  0 

xx  yy  zz 


(1) 


(ii)  Bernoulli's  equation  on  the  free  surface 


£ 

2 


2  2  2 

[  (4>x~U)  +4>y+())z  ]  +  Pgy  +  p  =  const. 


(iii)  on  the  body  surface,  z  =  f(x,y) 


f 

x 


(vu) 


+  f 


y  y 


(2) 


(3) 


(iv)  on  the  free  surface,  y  =  n(x,z) 


«>X-U)  -  ^y 


+  n 


=  o 


(A) 


Equations  (2),  (3),  and  (A)  are  exact  boundary  conditions.  In  order  to  derive 
the  solution  of  the  thin  ship  theory,  it  is  necessary  to  linearize  these  equations 
with  the  assumptions  that  the  derivatives  of  <f» ,  f,  and  n  are  small  and,  hence, 
second  degree  terms  are  negligible.  The  linearized  conditions  become 
(ii)  free  surface  condition 


P  +  k  <p  =  0 

xx  o  y 


(5) 
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(iii)  body  boundary  condition 

0  =  U  f 

z  x 

(iv)  the  free-surface  wave  height 


n 


f  4>x  (x.0,0) 


(6) 


(7) 


The  solution  of  Equation  (1)  with  the  boundary  conditions.  Equations  (5)  and  (6),  is 
given  by 


'Mx.y.z) 


U 

2tt 


G(x,y,z;x",y",0) 


f 

x 


(x',y")  dx'dy 


(8) 


where  G  is  the  Green  function  as  given  in  Reference  2  as 


G (x , y , z ; x  ,y  ,z')  = 


i  i  4k 

-  +  —  +  — 
r  r^  tt 


f  2  r 

I  sec  t>  dt*  I  dk  exp  [k(y+y  ) 
0  0 


cos  [k(x-x^)cosO ]  cos  (k(z-z^)sintV 


k  -  k  sec“0 
o 


J 


+  Auk  |  dO  sec  0 
o 


*  exp  [k^Cy+y^sec  G]  sin  [k^  (x-x  *)  secO  ] 


cos  [k  (z-z ') sinG • sec^G) 
o 


(9) 
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where  k  =  g/U 

r  =  [ (x-x')2+(y-y")2+(z-z")2]1/2 
r1  =  [  (x-x')"+(y+y')2+(z-z>)2]1//2 
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Here,  (x,y,z)  is  the  point  where  the  potential  is  sought  and  (x' ,y' ,z')  is  a  source 
point.  With  substitution  of  Equation  (8)  into  Equation  (7)  and  integrating  by  parts, 
the  free-surface  wave  height  becomes 


L 


„2 

(  ,,  r  1 

>  ' 

U 

Tl  2rrg 

|  JJ  G^x,0,0;x,y",0^fx"x'‘dx>dy  “  JGfx"dy" 

L  ; 

(10) 


The  fact  that  '  =  ~  Has  been  applied  in  the  derivation  of  Equation  (10). 


NUMERICAL  COMPUTATION 

In  order  to  simplify  the  numerical  computation.  Equation  (9)  is  expressed  in 
different  form' 


G  (x  ,  0,  0;  x  ”,  y  ' ,  0)  =  I  +  J 


4k 


O  o 

J"  sec2,3  HoJ* 


( i+kQy 'sec20)  e  ^  dS. 


(£+k  y'sec^S)2  +  [k  (x-x ') sec20  ] 2 
o  o 


-  8H(x 


V)  k°  J 


sec  J  e 


k  y  sec 
o 


2, 


sin  [k  (x-x")sec0)  d0  (11) 
o 


where  H(x)  is  the  Heaviside  function,  defined  as  1  for  x  >  0,  0  for  x  <  0,  and  1/2 
for  x  =  0.  By  substitution  of  Equation  (11)  into  Equation  (10),  the  free-surface 
wave  height  is  expressed  as 


4 


(I+J)f 


:'X'd  x>'-J 


(I+J)f  -dy " 
x 


The  first  term  of  Equation  (11)  can  be  rewritten  as  given  in  Reference  4 


2 

•/  d-i 


(?cos~0+k  y  ”)  e  d? 
o 


“*q  Jq  (Hcos^d+k^y  +  [kQ  (x-x')cosO 


where 


Q(a,b,0)  d0 


O 

i.b.O)  =J* 


(£cos^0+a)  e  V  dt 
(£cos"0+a)  ~  cos^f 


with  a  =  k^y '  and  b  =  kQ(x-x”).  The  general  behavior  of  the  integrand  of  Q(a,b,0) 
is  explained  in  detail  in  Reference  4.  The  method  of  numerical  evaluation  of 
Q(a,b,0)  in  this  report  is  different  from  that  of  Reference  4.  In  this  represen¬ 
tation  the  integrand  is  differentiated  with  respect  to  ?  and  set  to  zero  in  order  to 

■i 

find  extreme  values  with  <cos~0  +  a  =  v,  the  maximum  or  minimum  values  of  the 
integrand  occur  in  the  following  solution 


3  2,  2,2  2,  ,2  4n 

v  +  cos  0*v  +  b  cos  u*v  -  b  cos  0=0 


There  are  three  complex  solutions  of  Equation  (15)  which  may  be  designated  as  v^, 
v^,  and  where  the  maximum  or  minimum  values  occur  in  the  integrand  of  0. 
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Only  Che  real  solutions  are  used  in  the  numerical  evaluation,  and,  in  most  cases, 

there  will  be  one  or  two  real  solutions.  if  there  are  two  real  solutions,  these 

2 

two  solutions  are  separated  from  the  point  Z  =  -a/cos^O.  At  this  point  the 
integrand  of  Q  becomes  zero.  This  case  is  represented  graphically  in  Figure  1. 


Figure  1  -  Integral  Regions  of  Q 


For  example,  if  there  are  two  solutions,  the  integral  region  is  subdivided  into  four 

subregions:  A.  (0  <  Z  <  H  ) ,  A«  (£  <  t  <  i  ),  A  {i  <  i  <  ),  and  A,  (C,  <  ■.  <  °°) . 

1  —  a  2.  a  —  o  3o—  d  4  o 

He  re,  i  and  l.  are  determined  as  follows 
a  b 


where  K  is  a  constant  which  depends  upon  a,  b,  and  0.  The  specific  value  of  the 
constant  K  which  is  employed  in  any  particular  computation  falls  between  1  and  10 


6 


and  is  found  by  examining  the  numerical  stability  of  the  algorithm.  Each  hull  form 
and  operating  condition  may  require  the  setting  of  different  integration  limits. 

The  process  of  selection  of  K  is  straightforward  but  is  not  automatic,  and  requires 
the  attention  of  the  investigator.  When  a  and  b  are  small,  the  integrand  of  Q 
increases  slowly  between  0  and  f,^,  and  then  decreases  rapidly  between  and  , 
Between  i.  and  the  integrand  increases  rapidly  and  decreases  slowly  when  >.  is 

larger  than  Therefore,  the  discretized  distances  for  numerical  integration  are 

different  in  those  subregions.  In  subregion  A ^  and  A^,  for  example,  the  integration 
step  is  smaller  than  that  of  A^  and  A^ .  This  approach  is  employed  in  order  to 
reduce  computational  time. 

In  each  region,  the  integrand  of  Q  is  computed  numerically.  As  a  special  case, 
when  v  =  tt/2.  Equation  (14)  becomes 


Q 


1 

a 


(191 


In  order  to  show  the  behavior  of  the  Q-function,  the  integrated  values  at  each 
region  are  given  in  Table  1  for  a  =  -0.12  and  b  =  1.0.  The  first  column  indicates 
the  O-variable  and  the  following  three  column  pairs  show  the  integrated  values  of  Q 
at  three  different  divisions:  A^,  A0 ,  and  .  The  last  column  is  the  total  Q-values 
at  each  given  0.  The  integral  of  the  last  column  with  respect  to  0  becomes  the 
result  of  I-integral  The  integrated  values  of  Q  are  small  positives  when  9  is  about 
70  deg  and  less  in  this  case.  For  higher  " ,  Q  becomes  large  and  negative.  The 
I-integral  is  always  negative. 

In  the  computational  procedure,  the  I-integral  is  substituted  into  Equation 
(12),  and  the  first  and  third  terms  are  computed  numerically. 

The  second  term  of  Equation  (11)  can  be  rewritten  with  change  of  variable, 
sec3  =  t  as  follows 


J  = 


-8k  H(x-x^) 
o 


1 


.  2 
t 


sin  [k  (x-x')t]  dt 


(20) 
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TABLE  1  -  NUMERICAL  EXAMPLE  OF  Q(-0.12, 


If  f  -  is  expressed  in  Equation  (12)  as 


f  >  -  =  A..  +  B  +  y"(C,.+D  x') 

xx  ij  ij  ij  ij 


(21) 


the  second  term  of  Equation  (12)  can  be  integrated  piecewise  over  the  x  y  plane 
and  becomes 


M-l  N-l 


fk-x-  dx'dy' =  _8  £  JL  H(x_°  “H  lAij+Bijx'+y'(cij+Dijx')] 
*  »  1_1  •_!  *0 


i-1  J-l 


+  ^(vDijy')-rf(vDijx')-rt£Dij 


(22) 


where 


,  ‘  J2 

«>  k  v  t 
o 
e 


—  cos  [k  (x-x')t]  dt,  p  =  2,4 

t  p}f^i 


(23) 


and 


>,s  “  j 


,  >  2 
w  k  y  t 
o 
e 


-  sin  [k  (x-x‘)t]  dt,  p  =  3,5 

tp 


(24) 


With  f  *(-L/2,y')  =  E.  +  F.y'",  the  fourth  term  of  Equation  (12)  can  be  integrated 
X  J  J 


as 


N-l  „ 


k- dy' =  _8  2* 

J  j=l 


-Jl,s(Ej+Fjy')+^  J3,s 


(25) 


L 

x  =-  2 


9 


In  Equation  (25),  the  Heaviside  function  is  satisfied  by  requiring  .1  ^  and 
to  be  zero  at  x”  =  x. 

While  J  and  J  include  liarmonic  functions,  the  numerical  evaluation  is 
p,c  p,s 

performed  in  each  cycle  of  the  periodic  function 


k  (x-x')  t  = 
o 


0,1,2... 


(26) 


When  n  =  0,  t  =  n / 2 k  (x-x’)  and  l  must  be  larger  than  1.  In  the  region  between 
o  o  u,y.y 

1  and  t  ,  the  substitution  of  variables  u“  =  t*"  -  l  applied  to  Equation  (23)  yields 


J 

P.c 


koy'(l+u2) 

p+1 

(1+0  “ 


cos  [k  (x-x’) 
o 


du 


(27) 


Equation  (27)  is  more  convenient  than  Equation  (23)  in  numerical  computations  for 

n  <  1,  but  for  n  >  1,  Equation  (23)  is  used  in  the  numerical  evaluation. 

When  k^y ’  is  small,  there  are  many  harmonic  oscillations  before  the  exponential 

function  decays  for  large  t.  In  this  case  the  rational-exponential  term  in 

Equation  (23)  is  replaced  with  a  simple  linear  function  so  that  the  integral  can  he 

done  easily.  With  a  =  k^y',  b  =  k  (x-x'),  tj  =  (2n+l/2)”  and  t.,  =  (2n+2+l/2 ) " , 

J  becomes 

P.c 


I 

I 


2  2 
at" 
e 


t . 


J  «1> 


B. 

i 


cos  bt  dt  =  |  (B,t+B9)  cos  bt  dt  =  -j—  (t,-!,) 

C1 


(28) 


With  the  same  analogy  and  in  this  case,  t^  =  2n''  and  t -  (2n+2)  ',  J  ;  is  given  by 


I 


C2  2  }2 
at  r 

sin  bt  dt  =  (Bj  t+B0)  sin  bt  dt  =  U 


B, 


(29) 
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In  each  interval  where  the  integrand  becomes  zero,  is  computed,  and  with  Equa¬ 
tions  (28)  and  (29),  J  and  J  are  approximated,  respectively.  Values  of  J 

p,cp,s  P,C 

and  J  for  a  =  -0.12  and  b  =  1.0  are  computed  and  given  in  Table  2.  The  first  row 
P » s 

shows  the  results  of  Equation  (27)  with  t  =  1.5708.  The  second  row  shows  the 

o 

numerical  integral  of  Equations  (23)  and  (24)  with  1.5708  <  t  <  7.85398.  As  shown 
in  Table  2,  the  integrated  results  of  Equation  (27)  are  dominant  in  the  whole  region 
of  computation. 

In  order  to  validate  the  numerical  results  of  Equation  (12),  two  hull  forms 
have  been  selected:  a  Wtgley  hull  and  an  SL7  container  ship.  The  Wigley  hull  is  a 
mathematical  hull  form  which  is  represented  as 

f(x,y)  =  |  ^1.0-  -j^l.O-  (30) 

where  L  =  20  ft,  B  =  2  ft,  and  H  =  1.25.  The  block  coefficient  of  this  hull  is 
0.444.  The  SL7  is  a  container  ship  with  a  large  bulbous  bow;  its  sectional  hull 
offsets  are  given  in  Table  3.  Its  length  is  880.5  ft;  width,  105.8  ft;  and  draft, 
34,1  ft.  The  block  coefficient  of  the  SL7  is  0.53. 

Figure  2  shows  the  wave  profiles  of  the  Wigley  hull  for  F^  =  0.266  and  0.452. 
The  wave  profiles  are  computed  at  the  centerplane  where  z  “  0,  and  the  experimental 
results  are  measured  at  the  hull  surface.  For  F  =  0.266,  the  numerical  results 

5  n 

agree  fairly  well  with  the  experiment.  For  F  =  0.452,  there  is  some  discrepancy 
near  the  bow  and  the  region  aft  of  the  midship  section.  As  shown  in  Reference  1, 
at  higher  Froude  numbers,  the  numerical  results  of  various  analytical  methods  do  not 
agree  well  with  the  experiment  in  the  bow  region.  This  indicates  that  the  effect 
of  trim  and  sinkage  might  be  important  when  the  speed  is  high. 

Figure  3  shows  the  computed  wave  profiles  for  the  SL7.  The  agreement  between 
the  numerical  results  and  experiment^  is  not  satisfactory  for  this  hull.  The 
large  bulbous  bow  might  be  the  cause  of  the  discrepancies.  The  thin  ship  theory 
assumes  that  hull  function  is  closed  at  the  bow  and  stern,  and  the  assumption  is 
violated  for  this  hull. 
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TABLE  2  -  NUMERICAL  EXAMPLE  OF  J  AND  J 

p,c  p,s 

k  y'  =  -0.12  AND  k  (x-x')  =  1.0 
o  o 


Range  of 
Variable  u  or  t 

m 

- 1 

2,c 

J. 

3,s 

m 

Jr 

J>,S 

0-1.21136 
Equation  (23) 

0.67517 

0.27498 

0.52425 

- 1 

0.23789 

0.42996 

1.5708-7.85398 
Equation  (20) 

0.16246 

-0.04389 

0.04810 

-0.00867 

0.01476 

Total 

0.83756 

0.23109 

0.57235 

0.22922 

0.44472 
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TABLE  3  -  OFFSET  OF  SL7  (in  feet) 


LONGITUDINAL  COORDINATE,  2x/L 


Wave  Profiles  of  SL7  for  F  =0.1  and 


CONCLUSIONS 

The  thin  ship  theory  has  been  applied  for  the  computation  of  the  free-surface 
wave  height  produced  by  a  ship  in  steady  forward  motion.  For  a  Wigley  hull,  the 
agreement  between  the  computation  and  experiment  is  fairly  good.  However,  the 
numerical  results  for  an  SL7  show  large  discrepancies  from  the  experiment.  It 
seems  that  the  inclusion  of  the  effect  of  a  bulbous  bow  and  transom  stern  to  the 
present  method  is  necessary.  Further,  the  effects  of  trim  and  sinkage  should  be 
included  in  future  developments. 
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